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ABSTRACT 

We present detailed predictions for the redshifted 21 cm signal from the epoch of reionization. 
These predictions are obtained from radiative transfer calculations on the results of large scale 
(100//i Mpc), high dynamic range, cosmological simulations. We consider several scenarios 
for the reionization history, of both early and extended reionization. From the simulations 
we construct and analyze a range of observational characteristics, from the global signal, via 
detailed images and spectra, to statistical representations of rms fluctuations, angular power 
£^ ■ spectra, and probability distribution functions to characterize the non-Gaussianity of the 2 1 cm 

\q \ signal. We find that the different reionization scenarios produce quite similar observational 

. signatures, mostly differing in the redshifts of 50% reionization, and of final overlap. All sce- 

narios show a gradual transition in the global signatures of mean signal and rms fluctuations, 
which would make these more difficult to observe. Individual features such as deep gaps and 
bright peaks are substantially different from the mean and mapping these with several arcmin- 
utes and 100s of kHz resolution would provide a direct measurement of the underlying density 
field and the geometry of the cosmological HII regions, although significantly modified by pe- 
culiar velocity distortions. The presence of late emission peaks suggest these to be a useful 
target for observations. The power spectra during reionization are strongly boosted compared 
to the underlying density fluctuations. The strongest statistical signal is found around the time 
of 50% reionization and displays a clear maximum at an angular scale of i ~ 3000-5000. We 
rS \ find the distribution function of emission features to be strongly non-Gaussian, with an order 

^ . of magnitude higher probability of bright emission features. These results suggest that obser- 

vationally it may be easier to find individual bright features than deriving the power spectra, 
which in its turn is easier than observing individual images. 

Key words: cosmology: theory — diffuse radiation — intergalactic medium — large-scale 
structure of universe — galaxies: formation — radio lines: galaxies 
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1 INTRODUCTION 

Observations of the redshifted 21-cm line of hydrogen are currently 
emerging as the most promising approach for the direct detection 
of the epoch of reionization, and possibly even of the preceding 
period, the Cosmic Dark Ages. Several large radio interferome- 
ter arrays are currently either operational, under construction or 
being planned that have the potential to detect this redshifted 21- 
cm emission. These projects include PAST 1 (already operating in 



north-western China), GMRT 2 (operating in India), LOFAR 3 (un- 
der construction in the Netherlands), MWA 4 (to be located in West- 
ern Australia), and SKA 5 . 

The observations will be complicated due to the combined ef- 
fects of (in order of increasing distance from us) possible radio fre- 
quency interference from terrestrial emitters (in particular in the 
87-108 MHz FM band), ionospheric fluctuations, the galactic fore- 
ground, and unresolved intergalactic radio sources, all of which are 
much stronger than the expected redshifted 21-cm signal. There- 
fore it is crucial to understand the characteristics of the reionization 
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21 -cm signal in some detail. This can help both in the planning of 
the experiments, so that they are optimized for the expected signal, 
and once data is available, can direct us in separating the strong 
foregrounds and interpreting the signal correctly. Predicting the 
unique signatures of reionization is quite hard. The ACDM frame- 
work for cosmological structure formation is now well-established, 
the fundamental cosmological parameters are increasingly better 
constrained, and the basic physical processes during reionization 
are fairly well-understood. However, we will be probing redshifts 
for which little or no observational data is currently available, and 
many of the relevant parameters are at best poorly constrained. For 
example, a number of choices remain open regarding the nature 
of the ionizing sources during reionization, their numbers, photon 
production efficiencies and spectra. 

Another complication is the large dynamic range required 
from any simulation which aims to predict the redshifted 21 -cm 
emission. Within the hierarchical structure formation paradigm the 
ionizing photon emissivity during reionization is dominated by nu- 
merous dwarf-size galaxies, rather than by larger ones, which at 
high redshifts are too rare to make an appreciable contribution. On 
the other hand, the strong source clustering at high redshifts means 
that ionized bubbles quickly overlap locally and each H II region 
typically contains a large number of sources. This leads to the for- 
mation of large ionized regions, of size tens of comoving Mpc or 
more, which requires simulation volumes of size ~ 100 Mpc for 
proper treatment. Cosmological simulations resolving dwarf galax- 
ies in such large volumes are only now becoming possible. Even 
more difficult and computationally-intensive is the transfer of ion- 
izing radiation from the tens of thousands up to millions of indi- 
vidual galaxies found in such a large volume. A further challenge 
is the absorption of photons by collapsed objects, which happens 
at even smaller scales. This however, only becomes important after 
or around final overlap, when the opacity of the IGM to ionizing 
photons has dropped to a few percent and these (proto-)galaxies 
(e.g. Lyman Limit Systems) become the main source of opacity for 
ionizing photons. 

To date there have been only a few cosmological simula- 
tions which studied the redshifted 21 -cm. Most of these resolved 
this difficulty by considering small enough regions (thus achiev- 
ing the high resolution required) either during the Cosmic Dark 
Ages (Shapiro et al. 2005; Kuhlen et al. 2006) or during reioniza- 
tion (Carilli et al. 2002; Ciardi & Madau 2003; Gnedin & Shaver 
2004; Furlanetto et al. 2004; Valdes et al. 2006). An alternative ap- 
proach is to simulate large regions at low resolution and thus with- 
out resolving individual small sources, and to include their effect in 
an averaged, approximate way (Kohler et al. 2005). 

Recently we presented the first N-body and radiative trans- 
fer simulations which considered a sufficiently large volume, and 
at the same time resolved all dwarf galaxies in that volume and ac- 
counted for their individual contributions to reionization (Iliev et al. 
2006, hereafter Paper I). We first performed an N-body simulation 
with 4.3 billion particles in a (100 h~ x Mpc) 3 volume. This pro- 
vided us with all halos with masses above 2.5 x 1O 9 M and the 
cosmological evolution of the density field. We then imported these 
results into a new fast and precise radiative transfer code called C 2 - 
Ray (Mellema et al. 2006). In Paper I we presented our simulation 
methodology, as well as the results for the large-scale geometry 
(often also referred to as topology) of reionization, i.e. the size and 
number distributions of the H II regions in space, their evolution in 
time, and the power spectra of the resulting neutral and ionized den- 
sity fields. We showed that the fluctuations in both the neutral and 
the ionized density are strongly boosted due to the patchiness of 



reionization. We also demonstrated that the probability distribution 
function (PDF) of the ionized fraction and the ionized density are 
generally strongly non-Gaussian at all scales and quantified their 
level of departure from non-Gaussianity for the first time. We also 
derived the PDF and the mean optical depths to Ly-a photons (re- 
lated to the Gunn-Peterson effect in spectra of distant sources). In 
this paper we present detailed predictions for the redshifted 21 -cm 
signal, including results from several new simulations in addition 
to the original one from Paper I. We focus on generic predictions 
of the observational properties from our simulations, rather than 
aiming our analysis to specific radio-interferometry arrays. 

The lay-out of this paper is as follows. In Sect. 2 we describe 
the basic procedures of extracting the 21 -cm signal from the data 
and the assumptions that go into that. In Sect. 3 we discuss our 
simulations and their basic parameters and features. In Sect. 4 we 
present the evolution of the mean 21 -cm signal and its implications. 
Next, in Sect. 5 we discuss how the reionization geometry would 
be seen at the redshifted 21 -cm line and the dependence of this on 
the adopted observational beam shape. The first 21 -cm signals to be 
detected would quite possibly be individual, bright features. These 
are discussed in Sect. 6. Finally, the statistics of the 21 -cm emission 
signal is discussed in Sect. 7. Our conclusions are summarized in 
Sect. 8. 

Throughout this paper we use the concordance flat (Qk — 
0) ACDM cosmology with parameters (Q m , Ha, Q&, h, as, n) = 
(0.27, 0.73, 0.044, 0.7, 0.9, 1) (based on the first year WMAP data 
Spergel et al. 2003), where ft m , Qa, and Qb are the total matter, 
vacuum, and baryonic densities in units of the critical density, ag 
is the standard deviation of linear density fluctuations at present on 
the scale of 8/i _1 Mpc, and n is the index of the primordial power 
spectrum of density fluctuations. 



2 THE REDSHIFTED 21-CM SIGNAL 

The 21 -cm radio line, in either emission or absorption is due to 
a spin-flip transition of neutral hydrogen atoms. This transition 
quickly enters into equilibrium with the CMB photons and hence 
it needs to be decoupled from the CMB by some mechanism in 
order to become observable. This could occur either through col- 
lisions with other hydrogen atoms and free electrons (Purcell & 
Field 1956; Field 1959; Allison & Dalgarno 1969; Zygelman 2005) 
or through Ly-a photon pumping (Wouthuysen 1952; Field 1959; 
Hirata 2006; Chuzhoy & Shapiro 2005). 

The differential brightness temperature with respect to the 
CMB of the redshifted 21 -cm emission is determined by the den- 
sity of neutral hydrogen, pm, and its spin temperature, T s (see e.g. 
Morales & Hewitt 2004; Shapiro et al. 2005, for detailed discus- 
sions). It is given by 



err Ts - Tcmb n — T\ 

STb= l + z (1 " 6 } 



(1) 



where z is the cosmological redshift, Tcmb is the temperature of 
the CMB radiation at redshift z, and the optical depth r is (e.g. Iliev 
et al. 2002): 



r(z) = 



3AgAi T*nifi(z) _ 0.28 (l + z 



32ttT s H(z) 



0.28 (1 + z\ 3 / 2 ^ rx 



(2) 



where Ao = 21.16 cm is the rest-frame wavelength of the line, 
Aio = 2.85 x 10" 15 s _1 is the Einstein A-coefficient, T* = 
0.068 K corresponds to the energy difference between the two lev- 
els, 1 + 6 = phi/ (ph) is the mean density of neutral hydrogen in 
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units of the mean density of hydrogen at redshift z, and H(z) is the 
redshift-dependent Hubble constant, 

H(z) = HqE(z) 

= H [Q m (l + zf + Q k (l + zf + Ov] 1/2 (3) 

« i/o^ /2 (l + ^) 3/2 , 

where E(z) is the redshift dependent part, as defined here, Ho is 
the value at the present day, and the last expression is valid for 
z > 1. 

Assuming r « 1 (which is always correct for the redshifts 
of interest here, except for some lines of sight through mini-halos 
(Iliev et al. 2002), equation 1 becomes 

ST„ * (3.1mK)ft- l (1 + f ( Ts - TcMB > (l + J) (4) 
E(z) Ts 

« (27 mK )(i^) 1/2 (Zi^)(l + 5 ), 

with h being ifo in units of 100 km s _1 Mpc -1 . 

To obtain an appreciable signal, the spin temperature Ts 
should differ significantly from Tomb- The gas kinetic tempera- 
ture itself is expected to be above the CMB temperature due to 
heating by shocks, X-rays, and Ly-a photons, thus the redshifted 
21 -cm is generally in emission. Decoupling through collisions is 
efficient only at very high-z and in significantly overdense and 
heated regions, mostly inside collapsed halos (Iliev et al. 2002, 
2003; Shapiro et al. 2005). Once the first stars turn on, the Lya 
photons produced by them can easily pump the population to the 
gas kinetic temperature (e.g., Ciardi & Madau 2003). Thus, during 
reionization one can assume that globally Ts ~ T gas ^> Tcmb- 
In this case the precise value of Ts becomes irrelevant, and the ob- 
servable differential brightness temperature is only a function of 
redshift and H I density. This is the assumption we adopt through- 
out this paper. One should note that during the earliest stages of 
reionization, when sources were few and could only affect their lo- 
cal environments this assumption might not hold (Barkana & Loeb 
2005b; Chuzhoy & Shapiro 2005), but such effects are still difficult 
to implement in simulations. 



3 THE SIMULATIONS 

The details of our methodology were described in Paper I. Here 
we only summarize the relevant parameters and briefly discuss 
the main features of each of the simulations used in the cur- 
rent work. Our base N-body simulation was performed using the 
particle-mesh code PMFAST (Merz et al. 2005) and has a total of 
1624 3 = 4.3 billion particles on a 3248 3 computational mesh. 
The computational volume is (100 /i -1 Mpc) 3 . We found all ha- 
los with masses above 2.5 x 1O 9 M (corresponding to 100 par- 
ticles or more) on a large number of density time- slices. We then 
imported these results into a new fast and precise radiative trans- 
fer code called C 2 -Ray which we have developed (Mellema et al. 
2006). The code has been tested in detail against available analyt- 
ical solutions (Mellema et al. 2006) and in comparison with other 
radiative transfer codes (Iliev et al. 2006). The radiative transfer 
runs have mesh resolutions of 203 3 and 406 3 , since the full N-body 
mesh, at 3248 3 , is still well beyond the current computational ca- 
pabilities. Our ionizing sources correspond to all the halos found in 
our volume at the full resolution of the underlying N-body simula- 
tion. No artificial grouping of sources has been employed, with the 
exception of combining sources which happen to be inside the same 




Figure 1. A cut through the simulation box of our high resolution, 
f2000_406 run, at z= 13. 89. The neutral gas density is shown in green. The 
orange overlay shows the ionized regions. Each side of the square corre- 
sponds to a comoving length of 100 Mpc. 



Table 1. Simulation parameters and global reionization history results 
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radiative transfer cell, something that affects only a small fraction 
of the sources. 

We present the results of five simulations, all of which use the 
same box size, density fields and halo catalogues, but each with 
different the sub-grid physics, as follows. Four of these simulations 
are run at the same 203 3 mesh resolution, but making different as- 
sumptions about the gas clumping at small scales and about the 
photon production efficiencies of the ionizing sources. The fifth 
simulation is one run at the higher mesh resolution of 406 3 , but 
otherwise adopting the same assumptions. All simulations and their 
basic parameters and features are summarized in Table 1 . 

Our first case, which we will call f2000 hereafter, is the simu- 
lation we presented in detail in Paper I. It assumes a source photon 
production efficiency (which is a combination of the total number 
of ionizing photons emitted by the stars per unit time, the star for- 
mation efficiency and the escape fraction) of f 1 — 2000 photons 
per halo atom, which corresponds to a top-heavy initial mass func- 



4 G. Mellema, et al. 



tion (IMF). It reaches overlap (defined as more than 99% ionization 
by mass) at a redshift of z — 11.3. The resulting optical depth for 
electron scattering is r es = 0.145, within the 2-cr limits of the new 
WMAP value, r es = 0.09 ± 0.03 (Page et al. 2006). Our second 
case is the same simulation as f2000 but at a higher mesh resolution 
of 406 3 (hereafter called f2000_406). In general the results closely 
match the f2000 run. The density field is resolved better and the H II 
regions have somewhat less spherical shapes than in the lower res- 
olution run (as could be expected). The global evolution is slightly 
slower due to the better-resolved density field, resulting in higher 
effective gas clumping. As an illustration we show in Fig. 1 a slice 
through the simulation volume of the density and ionization struc- 
tures at redshift z — 13.89, extracted from this high-resolution run. 

Our third simulation (labelled f250) adopts a lower photon 
production efficiency of 250 photons per halo atom, corresponding 
to either a slightly top-heavy IMF, or a Salpeter IMF combined with 
a bit higher escape fraction and star formation efficiency at these 
high redshifts. This simulation reaches overlap at redshift z = 9.3. 
The resulting optical depth for electron scattering is r es = 0.121, 
within the l-a limit of the new WMAP value. 

The fourth and fifth cases included in this paper, called f2000C 
and f250C, assume the same source efficiencies as f2000 and f250, 
but add the effect of sub-grid gas inhomogeneities, described by 
a mean volume- averaged clumping factor, C su b g rid = (n 2 ) /(n) 2 , 
given by 



C snbetid (z) = 27.466e-°' 114z+0 ' 001328 z2 . 



(5) 



This fit to the small-scale clumping factor is an improved version 
of the one we presented in Iliev et al. (2005). To derive it we used 
another PMFAST simulation, with the same computational mesh, 
3248 3 , and number of particles, 1624 3 , but a much smaller com- 
putational volume, (3.5 h -1 Mpc) 3 , and thus much higher resolu- 
tion. These parameters correspond to particle mass of 1O 3 M and 
minimum resolved halo mass of 1O 5 M . This box size was cho- 
sen so as to resolve the scales most relevant to the clumping - on 
smaller scales the gas would be Jeans smoothed, while on larger 
scales the density fluctuations are already present in our computa- 
tional density fields and should not be included again. The expres- 
sion in equation (5) excludes the matter residing inside collapsed 
halos since these contribute to the recombination rate differently 
from the unshielded IGM. The minihalos are self- shielded, which 
results in their lower contribution to the total number of recombi- 
nations than one would infer from simple gas clumping argument 
(Shapiro et al. 2004; Iliev et al. 2005), while the larger halos are 
ionizing sources and their recombinations are implicitly included 
in the photon production efficiency f 1 through their escape frac- 
tion. The f2000C case reaches overlap at a redshift of 10.15, and 
the case f250C at 8.2. The resulting integral optical depths for elec- 
tron scattering are r es = 0.135 and 0.107. 

The value for the electron scattering optical depth reported 
from the 3-year WMAP data, 0.09 ± 0.03 (Page et al. 2006), is 
substantially below the first year value (0.17 ± 0.04, Kogut et al. 
2003). Our simulation results fall roughly within the overlap region 
between these two results and within 2-cr from the new result. Mod- 
erate changes in the simulation parameters, e.g. assuming lower lu- 
minosity sources and/or accounting better for the small-scale gas 
clumping as in f2000C and f250C easily extends the reionization 
process until significantly later times. However, the value of the 
optical depth should not be considered separately from the values 
of the underlying cosmological parameters (Spergel et al. 2006). 
Alvarez et al. (2006) showed that to first approximation simula- 
tions employing the first year WMAP cosmological parameters and 




Figure 2. Simulation size and output times: (top) the frequency of the 21- 
cm signal (the crosses indicate the redshifts at which we have model data), 
(middle) the angular size on the sky of our computational volume, and (bot- 
tom) the frequency width of our computational volume, all plotted vs. red- 
shift z. 



achieving a r of 0.17, would for the 3-year WMAP cosmological 
parametes achieve a r of 0.1. Based on these results, we expect that 
our current results would all be consistent with the 3 -year WMAP 
data. This could be checked with future simulations. 

In this paper we will use results from simulations, f2000, 
f2000C, f250 and f250C, and results from f2000_406 only as com- 
parison to f2000. A more detailed comparison of their different 
reionization histories and implications will be the subject of a fu- 
ture paper. 

For reference we show in Fig. 2 the frequency of the redshifted 
21 -cm emission, v z (top), the angular size on the sky, A#box (mid- 
dle), and its total bandwidth, A z/ box (bottom), versus redshift, for 
our 100 /i -1 Mpc 3 comoving computational volume. These quan- 
tities are given by 



A# bo 



Az/ bo 



i/o/(l + *) 

^box 

(l + z)D A (z) 

bo 



(6) 
(7) 

(8) 



c(l + zY ' 

where Da(z) is the angular diameter distance, Lbox = 
100 /h Mpc is the comoving size of our box, c is the speed of light 
and uq — 1.420 GHz is the rest-frame frequency of the line. We 
see that our computational box corresponds to almost one degree 
on the sky and a range of 5.5 to 7 MHz in frequency. The intrin- 
sic resolution of our data is approximately 14" (7") in the spatial 
direction, and 30 KHz (15 kHz) in frequency for our 203 3 (406 3 ) 
meshes. None of the currently planned experiments will achieve a 
similar spatial resolution. The estimates for both PAST and LOFAR 
are that they will have synthesized beams of ~ 3'. Their intrinsic 
frequency resolution will actually be somewhat higher than in our 
simulations, e.g. around 1 kHz for LOFAR. However, to obtain suf- 
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Figure 3. Evolution of (top) the mean mass ionized fraction, x m , and (bot- 
tom) the mean differential brightness temperature, 5Tb, both as a functions 
of frequency for f2000 (solid, blue), f2000C (short-dashed, red), f250 (long- 
dashed, green), f250C (dotted, magenta). 



ficient flux sensitivity, the data will be integrated over significantly 
larger band widths, e.g. around 200 kHz for LOFAR, which approx- 
imately corresponds to 7 (14) of our cells at 203 3 (406 3 ) resolution. 



4 THE EVOLUTION OF THE MEAN BACKGROUND 

The simplest global measure of the progress of reionization is the 
evolution of the globally-integrated mean 21 -cm signal. In Fig- 
ure 3 we show how the mass-weighted ionized fraction x m and 
the mean differential brightness temperature, 5Tb evolve against 
observed frequency, zy obs for f2000, f2000C, f250 and f250C. The 
first thing to note is that the reionization is significantly delayed, 
by Az ~ 1 compared to f2000, when the small-scale clumping of 
the gas is included, and even more so, by A z ~ 2, when the ion- 
izing sources are less efficient photon producers. Accordingly, the 
mean signal disappears and becomes undetectable at frequencies 
above ~ 110 MHz in f2000, but does so at much higher frequen- 
cies, - 120 MHz (~ 135 MHz, - 145 MHz) in f2000C (f250, 
f250C). At face value this seems to imply that our f2000 reioniza- 
tion history will be impossible to observe in places with high FM 
interference, as is the case with e.g. LOFAR and GMRT, which op- 
erate above 115 MHz, but as we will show below, this is probably 
a too simplistic conclusion. 

Note that in the real universe the signal will never go com- 
pletely to zero since there will be neutral hydrogen contained in 
collapsed objects such as damped Lya systems and Lyman limit 
systems. However, at the redshifts we are considering these sys- 
tems will never contain more than a few tenths of a percent of the 
mass in our computational volume, and their combined signal will 
be too weak to be observable at 21 cm. 

It has been proposed that one could possibly detect the tran- 
sition from fully-neutral gas to almost fully-ionized as a "global 



reionization step" over the whole sky (Shaver et al. 1999). To be de- 
tectable against the smoothly varying foregrounds, such a transition 
should occur fairly quickly, resulting in a sharp drop in the global 
21 -cm signal in frequency space when observed with a single dish 
radio telescope. All of our models show a rather gradual transition, 
with the mean signal decreasing by ~ 20 mK over ~ 20 MHz. De- 
tecting such a transition is in principle well within the capabilities 
of even a 10 m radio dish (Shaver et al. 1999). However, a gradual 
change over 20 to 30 MHz would be difficult to disentangle from 
the strong foregrounds. The practicality of such observation would 
thus depend strongly on the (still highly-uncertain) detailed prop- 
erties of the spectral index variations of the foregrounds and how 
well one can model and subtract them. 



5 REIONIZATION GEOMETRY SEEN IN REDSHIFTED 
21-CM EMISSION 

5.1 Global evolution 

We start with an overview of the progress of reionization in our sim- 
ulations f2000 and f250, as seen in the 21 -cm emission line. We 
present a visualization of the evolution of the differential bright- 
ness temperature in Fig. 4. We constructed these slices by linearly 
interpolating the evolution of the whole simulation box between 
our discrete time outputs, wrapping the box periodically after each 
complete crossing of the volume, and taking slices through the re- 
sulting long box. The vertical axis corresponds to the complete ex- 
tend of our simulation volume of 100 Mpc comoving, while 
the horizontal axis shows the redshift at that position. Since red- 
shift corresponds to frequency, these slices are equivalent to slices 
through an idealized observational image-frequency volume. We 
took the slices at an angle to the simulation cube axes, to avoid ar- 
tificial periodicity effects caused by repeatedly passing through the 
same structures at later times. We also include the effect of pecu- 
liar bulk velocities from our N-body simulations, which result in 
redshift distortions along the line-of-sight (see Sect. 6.2 for more 
discussion on this last point). 

Both slices show all the basic evolution trends we discussed in 
Paper I. The reionization starts around z ~ 20, with a few isolated 
and highly clustered ionized regions. Most of the gas is still neutral 
and quite bright at 21 -cm emission, with some pixels reaching 5Tb 
of well over 100 mK. The Cosmic Web is already well-developed 
and shows in the images even at these fairly large scales. For f2000 
(Figure 4, top image) the first H II regions quickly merge with each 
other locally, forming significantly larger ones, of size ~ 10 Mpc 
comoving, already by z ~ 14, while at the same time many new 
ionizing sources emerge and start expanding their own local H II 
regions. By redshifts z ~ 12.5 many of the larger H II regions 
merge together to form a couple, and later on just a single, very 
large and topologically-connected (but highly-nonspherical) vol- 
ume of > 10 4 (Mpc) 3 . Essentially all remaining ionized regions 
merge with this large one around z ~ 11.5 and the last pockets 
of neutral gas (mostly in the large voids, which have no local ion- 
izing sources and are fairly far away from any source) are ionized 
by z ^ 11. In the f250 case (Figure 4, bottom image) the evolu- 
tion is a bit slower due to the weaker sources, and the H II regions 
initially smaller, on average. Essentially the same main evolution- 
ary trends remain in force, with significant local overlap of H II 
regions starting at z ~ 12, and the wide-spread merging of these 
between redshift z = 11 and 11.5. These finally merge into one 
huge ionized region between z — 10 and 9.5, although significant 
local pockets of neutral gas remain until the end. 
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Figure 4. Position-redshift slices from the f2000 (top) and f250 (bottom) simulation. The spatial scale is in comoving units. These slices illustrate the large- 
scale geometry of reionization seen at 21 -cm and the significant local variations in reionization history. Observationally they correspond to slices through an 
image-frequency volume. 



While these features are generic and thus seen in all slices, 
there are a number of special features to note. Reionization is a 
highly inhomogeneous process, with large variation in its progress 
between different 2D slices of the same simulation, and even more 
between different lines-of- sight (LOS). For example, in the f2000 
slice some LOS contain only few 21 -cm features after z ~ 13, 
while others go through neutral pockets even at global overlap 
(z 11). In the f250 case similar LOS behaviour can be seen. 
This strong patchiness until late times will have important impli- 
cations for direct observations of high-z Ly-a sources, a point we 
will address in future work. 



5.2 Sky maps 

In what follows we approximate the synthesized beam of angular 
size AO with a compensated Gaussian filter, given by 

^coW = ^(l-^)exp(l-^), (9) 

with a FWHM equal to AO = 2a^/2(l - LambertW(e/2)) 
(LambertW(e/2) « 0.685). The compensated Gaussian has a 
shape in Fourier space whose real part, k 2 exp(— k 2 /2), approx- 
imates well the actual observational beam shape of a compact 
interferometer (often referred to as 'dirty beam'), being insensi- 
tive to large scale features. For this reason we use the compen- 
sated Gaussian in most of our analysis. Its average value is zero, 
with a Gaussian- type peak in the middle, surrounded by a nega- 
tive through. As a consequence, data convolved with this beam will 
have both positive and negative values. In some cases we also use 
two other beam shapes, a Gaussian with FWHM of AO and a top- 
hat function of width AO. These two beam shapes are simple and 
widely used in the literature. We use them in a few cases below, 
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log 10 (l) 

Figure 5. Angular power spectra of the compensated Gaussian (red/solid), 
Gaussian (short-dashed/blue), and tophat (long-dashed/green) beams all 
with FWHM of 3 arcmin and each normalized to maximum of 1. The 
sharp cut-off of the tophat beam produces high t oscillations, showing 
that it should not be used for Fourier analysis of results. The compensated 
Gaussian, which is the one most closely resembling real-life interferometer 
beams has less power at low ts than the Gaussian, reflecting the relative 
insensitivity of compact interferometers for large scale structures. 



and when we do so we explicitly specify it. For reference we show 
plots of the angular power spectra of all three beam types in Fig. 5. 
The frequency bandwidth integration is always done with a tophat 
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Figure 6. The geometry of reionization seen at redshifted 21-cm line at redshifts z = 16.08, 13.62 and 12.57 for f2000. Top row panels: a thick slice 
(corresponding to 0.2 MHz in frequency) showing H II regions (blue) superimposed on the density field (green). The ionized density increases from dark 
blue to white, the neutral density from dark green to white. The dark green and dark blue patches are those areas in the thick slice that contain a mixture of 
ionized and neutral material. Middle row panels: The same reionization topology as seen at 21-cm with 3 / FWHM compensated Gaussian beam and 0.2 MHz 
bandwidth. Bottom row panels: The same as the middle panels, but with a Gaussian beam of the same FWHM. 



function, which is generally a good approximation to the actual in- 
tegration used in observations. 

In Fig. 6 we show sample 'thick' slices of the density fields 
(green/white) with the ionized regions shown in blue, at three red- 
shifts, taken from f2000 A thick slice contains the average of sev- 
eral slices corresponding to a given band in frequency, to allow 
comparison to the images with finite bandwidth below. We show 
thick slices for redshifts z = 16.08, early in the evolution (mass- 
weighted ionized fraction x m — 0.026), at z = 13.62 when the 
volume is about half-ionized (x m = 0.52), and at at z = 12.57, 
when most of the volume is already ionized, but significant neu- 
tral gas patches remain (x m = 0.80). Below each of those slices 
we show the corresponding 21-cm emission maps, each integrated 



over a bandwidth Az/=200 kHz and convolved with a compensated 
Gaussian beam of FWHM 3'. The third row shows the images for 
a Gaussian beam of the same FWHM. For simplicity, when inte- 
grating over the bandwidth we neglected the redshift evolution and 
distortions over that bandwidth, which is justified for such small 
Aza As explained above, the compensated Gaussian beam pro- 
duces negative values since its average is zero. 

An important point to note is that while the minimum and 
maximum values change over time, the full range from the bright- 
est to the least bright pixel in this redshift interval is roughly inde- 
pendent of the redshift and remains around 30 mK for a Gaussian 
beam and around 40 mK for a compensated Gaussian. However, 
this is not the case very early, when the reionization has barely 
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Figure 7. (top) Sample 21-cm sky maps at redshifts z = 11.68 and 10.14 from f250 and (bottom) the same maps, but after subtracting a wide, 5 MHz band 
from it. Subtracting the signal from such a wide band should remove most of the foregrounds without affecting the underlying redshifted 21 cm signal. 



started and much later, close to overlap; in at those times the fluc- 
tuations are less pronounced. This happens because during these 
epochs the fluctuations are mostly dictated by the fluctuations of the 
underlying density field, rather than by the patchiness of reioniza- 
tion, which strongly dominates the fluctuations at the intermediate 
times. The images illustrate how well the density fluctuations and 
the H II regions would show on a realistic 21-cm sky map. In spite 
of the significant level of smoothing, all large-scale features remain 
clearly visible in the maps. The larger H II regions are apparent as 
holes in the 21-cm radiation, and the neutral density peaks corre- 
spond to strong emission peaks. Some smaller, isolated ionized re- 
gions, with sizes below the beam smoothing length disappear from 
the map, particularly for the Gaussian smoothing. However, both 
the density fluctuations in the neutral regions and the large-scale 
patchiness of reionization remain clearly visible. The compensated 
Gaussian has much less power on large scales than the Gaussian 
beam, but more power on smaller scales, and therefore shows the 
detailed features better. The weaker 21-cm signatures of the mostly- 
ionized structures inside the H II regions are washed out and would 
require significantly higher spatial resolution and sensitivity in or- 
der to be detected. 

Although obtaining images of 21-cm emission would be the 
ideal observational result it is not expected that such maps will be 
easily extracted in the near term, due to the sensitivity limits of the 
planned observations, and especially due to the strongly dominant 
foregrounds. In order to eliminate the foregrounds effectively one 
has to make use of their frequency characteristics. The issues of 



foregrounds is a complicated one, and goes beyond the scope of 
this paper. However, if one naively approximates the foregrounds 
as a slowly-varying power-law in frequency with smooth angular 
variations, we can simply subtract the mean signal over a given 
band Aza Here we take Av to be 5 MHz. Such an approach would 
average out any small random variations of the spectral index, as 
well as subtract the main, power-law-like component. 

Figure 7 shows two sample maps (using a compensated Gaus- 
sian beam; top images) and the corresponding maps with a 5 MHz 
band subtracted (bottom images). These correspond to the half- 
ionized point (z = 11.68, left panels) and the late phase (z = 
10.14, right panels) of our f250 simulation. The subtraction pro- 
cess leaves the structures in the map largely intact, demonstrating 
that such a subtraction procedure would not affect the signal sub- 
stantially, while at the same time eliminating much of the assumed 
foregrounds. How well such a procedure would work in practice 
would of course depend on how well-behaved the foregrounds are, 
and in reality more intricate methods surely are needed. 

Obtaining detailed 21-cm sky maps would give us a very rich 
set of data on both early structure formation and the progress of 
reionization. However, extracting these from the noisy data domi- 
nated by strong foregrounds will not be easy. Thus, in the following 
sections we first consider the detectability of individual large fea- 
tures and then statistical measures of the expected signal, both of 
which should be significantly easier to obtain from the observa- 
tional data, and so should be the first goal of the observations. 
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Figure 8. The brightest peak in the simulation box as a function of redshift, for (from left to right, simulations f2000, f2000C, f250, as labelled). 
Shown are the maximum pixel value of the differential brightness temperature, STb i7nax vs. redshift z for several beam-sizes and bandwidths, as follows, 
(A#beam 5 ^bw) — (6, 0.4) (blue, solid), (3,0.2) (red, short-dashed), (1,0.1) (green, long-dashed) and the full resolution of our simulation (i.e. correspond- 
ing to single simulation cells; magenta, dotted), where the beam sizes are in arcminutes and the bandwidths are in MHz. For reference we also show the mean 
differential brightness temperature over the whole box (thin, black, solid) and the same, but if the gas were fully-neutral (i.e. if no reionization occurred; thin, 
black, dotted). 



6 INDIVIDUAL FEATURES 
6.1 Rare, bright peaks 

The brightest emission peaks for a given beam and bandwidth are 
possibly the 21 -cm emission features easiest to detect. In Figure 8 
we show the maximum pixel value of the differential brightness 
temperature, ATb, ma x versus redshift z for beam-sizes and band- 
widths (A#beam, Avbw) = (6, 0.4) (roughly corresponding to the 
ones relevant to PAST), (3,0.2) ("LOFAR"), (1,0.1) ("SKA") and 
the full resolution of our simulation (corresponding to single cells), 
where the beam sizes are in arcminutes and the bandwidths are in 
MHz. In order to compare to the mean signal and the signal without 
reionization, we use a Gaussian beam here, since the compensated 
Gaussian has a zero mean by definition. 

The first thing to note is that the values and the shape of the 
curves do not depend significantly on the reionization history, ex- 
cept for the a shift in redshift, by about Az = 1.5 (2) between 
f2000 and f2000C (f2000 and f250). The magnitudes of the bright- 
est peaks is fairly high, at several tens of mK even in the beam- and 
bandwidth- smoothed cases, and well over 100 mK for full resolu- 
tion. The magnitudes generally drop with redshift, but only mildly 
and the curves stay largely flat over the whole redshift range. The 
higher resolution and sensitivity of SKA would provide a signifi- 
cant improvement, while the difference in the signal between the 
LOFAR and PAST scales is rather marginal. 

These brightest peaks generally correspond to high-density, 
still neutral regions, i.e. regions that are on the verge of form- 
ing galaxies, who in turn would start ionizing their surroundings. 
Hence, while the magnitude of the peaks do not vary much with 
redshift, the location in the box at which it is found does vary signif- 
icantly and the brightest peaks at one redshift become deep troughs 
at later times. For reference we also show the mean simulated 21- 
cm emission signal, as well as the mean signal if we assume the 
volume to be completely neutral (i.e. if reionization never hap- 
pened). The mean signal from the box is well below the maximum 
one at all redshifts, even for maximum smoothing, and it decreases 
much more steeply with redshift, thus the difference from the peaks 
grows significantly with time, from factors of 1.2-6 (depending on 
the smoothing) at z ~ 16, up to 35-320 at late times (z ~ 11.5). 
Notably, the peak signals are also generally higher than the mean 



signal for a neutral volume. This indicates that the peak signals al- 
ways come from overdense regions, even though, as we showed in 
Paper I, the reionization is clearly inside-out, i.e. the high-density 
peaks get ionized preferentially compared to the voids. For high 
levels of smoothing (i.e. large beam sizes and bandwidths) the peak 
signal drops below the corresponding fully-neutral gas signal since 
by that time our volume is more than half-ionized and hence any 
beam that samples much of the box would inevitably include sig- 
nificant ionized regions. This is a consequence of the high bias of 
the density peaks (and sources) at high redshift, which means that 
the maximum- signal peak inevitably is close to other density peaks 
that have been ionized already. 

Searching for the relatively rare, bright peaks could thus be 
an efficient way to detect the signatures of reionization, even in its 
late stages when the universe is already more than 99% ionized. 
This would make these stages observable to radio arrays which are 
strongly affected by interference in the FM-band, such as LOFAR 
and the GMRT, since they correspond to frequencies well above the 
FM band. The statistics of these rare events is discussed in Sect. 7.3. 



6.2 Line-of-sight spectra 

Next we present line-of-sight (LOS) 21 -cm emission spectra ob- 
tained from our simulations. These can in principle probe in de- 
tail the neutral structures at high-redshifts along that LOS, as well 
as detect the H II regions as deep troughs in the spectra. In order 
to obtain spectra over large frequency bands we use the technique 
discussed in Sect. 5.1 for deriving the evolution over the complete 
redshift and time intervals of our simulations. We present the full- 
resolution spectra (essentially just LOS cuts through the data shown 
in Fig. 4), as well as the corresponding beam- and bandwidth- 
smoothed spectra for our simulation f250, since its overlap is lat- 
est, thus it extends furthest into the most easily observable redshift 
range. 

The full-resolution (at resolution ~ 30 kHz, corresponding 
to one computational cell) spectra are shown in Figure 9. The red 
line shows the spectra including redshift distortions due to the bulk 
peculiar velocities, and the green line is the difference between the 
distorted and undistorted spectra. The redshift distortions introduce 
substantial changes in the spectrum, moving peaks to different fre- 
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Figure 9. Sample line-of-sight 21 -cm spectra obtained from our simulation data (simulation f250). The figure shows the redshift distorted spectra (red) and 
the difference between the distorted and the undistorted spectra (green). 



quencies (showing up in the difference as a deep dip next to a high 
peak), merging the emission of some gas parcels, and spreading 
others out over a larger range in frequency. The comparison shows 
that these are noticeable effects that should always be taken into 
account. These redshift distortions can in principle be used to de- 
rive information about the underlying cosmology as well as about 
reionization (Barkana & Loeb 2005a; McQuinn et al. 2005), but 
this falls outside the scope of the current paper. 

As was also noted above, at high resolution there are some 
very bright pixels, with 5Tb > 50 mK, and these are fairly com- 
mon during most of the evolution, and generally at least one is 
found on every LOS. However, they become increasingly rare after 
z ~ 12 (z^obs ~ 110 MHz), since by that time most of the high- 
density peaks have formed ionizing sources. Conversely, at about 
that time the H II regions become common and start to be seen 
along most LOS as deep troughs. The size of these troughs initially 
corresponds to the typical size of the ionized regions, ~ 10 Mpc co- 
moving, or about 0.5 MHz. As discussed in detail in Paper I around 
the time the volume is half-ionized these ionized regions quickly 
start merging with each other and start forming much larger ones, 
seen here are wider troughs, of several MHz each. Most of the sig- 
nal eventually disappears as the volume becomes ionized. However, 
as we already pointed out, this occurs at quite different frequen- 
cies along different LOS. Along some LOS there are essentially no 
neutral structures for v > 120 MHz, while along others significant 
neutral regions persist until ^ 135 MHz or later. 

These fairly sharply-defined H II regions suggest a local (on 
scales tens of Mpc) equivalent of the "Global step" discussed in 
Sect. 4. While the last proved rather gradual, on a more local level 



the steps are be quite steep, from 30-50 mK down to almost and 
back within a fraction of a MHz. 

In Figure 10 we show the effect of beam- and bandwidth- 
smoothing on the spectra. The beam we use is again the com- 
pensated Gaussian, which produces negative differential brightness 
values at the minima. Clearly the smoothing reduces the signal con- 
siderably when most of the material is still neutral {y < 100 MHz), 
and the fluctuations are entirely due to density fluctuations, at rela- 
tively small scales. Once the first H II regions begin to expand, they 
produce much stronger fluctuations. These show up in the spectra, 
although the smallest ones are still averaged out by the smoothing. 
The valleys corresponding to the H II regions are generally deeper 
than the neutral density peaks, so they should be easier to detect. 
At later stages the isolated peaks are suppressed in some cases, but 
often remain clearly visible, depending on their characteristic size 
and environment. For these LOS and smoothing the signal never 
varies more than a few tens of mK, but the main features of the 
unsmoothed data are preserved and should still be detectable. 



7 STATISTICAL SIGNALS 
7.1 Spatial rms fluctuations 

The prominent individual 21 -cm features we discussed in the pre- 
vious section, either rare, bright emission peaks or large individual 
H II regions would quite possibly be the first signatures of the reion- 
ization epoch to be seen. Once such a first detection is made the 
next aim would be to put more robust constraints on the progress 
and duration of reionization than the ones currently available. At 
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Figure 10. Sample line-of-sight 21 -cm spectra obtained from our simulation data (simulation f250). Shown are the full-resolution (red) and the beam- and 
frequency-smoothed spectra (blue). For the latter we used a compensated Gaussian beam with a FWHM of 3' and a bandwidth of 0.2 MHz. 



the next level of complexity, the data may be used to derive the sta- 
tistical measures of the fluctuations in the 21 -cm signal. This could 
be achieved by considering a number of observational fields, thus 
suppressing the noise and systematics in the foregrounds. 

The fluctuations of the 21 -cm EOR signal can come from sev- 
eral sources. Here we concentrate on the dominant effect, the large- 
scale patchiness of reionization. The fluctuations can be derived as 
a single number (rms), for any given scale (defined by the observa- 
tional beam and bandwidth), or in more detail in the form of angular 
power spectra. We consider both here, starting with the rms fluctu- 
ations. The rms fluctuations are the equivalent of the global step 
discussed in Sect. 4, but observed with an interferometer (which is 
insensitive to a mean global signal). 

In Fig. 11 we show how the fluctuations in the mean signal 
change with redshift/frequency (top panel) for f2000, f2000C, f250 
and f250C. These were calculated using a top hat beam of 3' and 
a bandwidth of 0.2 MHz. The top hat was used here for computa- 
tional convenience, using a different beam would change the num- 
bers slightly, but would not affect the basic trends which we want 
to illustrate here. Early-on the 21 -cm fluctuations closely follow 
the gas density fluctuations, since very little of the gas is ionized. 
Once some more substantial ionized patches start appearing, with 
sizes roughly of order of the beam- and bandwidth smoothing, 
the differential brightness temperature fluctuations quickly raise 
and reach a maximum of (5T b 2 ) 1/2 ~ 10 mK, after which they 
drop as more and more of the universe becomes ionized, reaching 
~ 1 mK at overlap. This behaviour is quite generic and occurs in 
all cases presented here. However, the frequency, ^ max (and red- 
shift) at which the maximum of the fluctuations is reached varies 



significantly among the runs, from ^ max ~ 97 MHz for f2000 up 
to z^max ~ 120 MHz for f250C. The peak value of the fluctua- 
tions decreases modestly for more extended reionization scenarios, 
from 11 mK for f2000, down to 7 mK for f250C. This is due to the 
fact that both clumping and lower photon efficiency move power 
to smaller scales, below the smoothing scale used here. In all three 
cases the peak fluctuations occur at the time when the universe is 
~ 50% ionized by mass. The peaks are roughly symmetric in all 
cases and significantly wider for the late-overlap scenario (f250) 
than for the early-overlap one (f2000). The increasing clumping 
C(z) with time, combined with an increasing number of sources 
leads to a broader peaks for the models with redshift dependent 
clumping (f2000C and f250C). 

As a reference we also show the STbdn 2 ) 1 ^ 2 /n), i.e. the fluc- 
tuations for reionization scenarios with the same STb(z) as the ac- 
tual simulations, but assuming that the process occurred homoge- 
neously rather than patchily. In this last case the fluctuations of the 
emission would have been solely due to the fluctuations in the un- 
derlying density field. Clearly the patchy reionization creates sub- 
stantial extra power over a wide range of frequencies. The differ- 
ence is especially large close to overlap, where the mean signal 
essentially disappears, while the actual patchy signal is ~ 1 mK 
even at overlap, when the mass neutral fraction is less than 1 per 
cent. The bottom panel underlines this point further by showing the 
rms fluctuations of the differential brightness temperature relative 
to its mean, (ST 2 ) 1 ^ 2 /STb, for all four cases, and the same for the 
density fluctuations, (n 2 ) 1 ^ 2 /n (where the last is the same for all 
cases, since they use the same underlying density field). We see 
that the relative differential brightness temperature fluctuations are 
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Figure 12. 2D angular power spectra of the differential brightness temperature fluctuations (solid lines) for f2000 (left) and f250 (right). For reference we also 
show the corresponding fluctuations if there were no reionization (i.e. due just to the fluctuations of the underlying gas density field assuming all the gas were 
neutral). 
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Figure 11. (top) Rms fluctuations of the differential brightness temperature, 
(8Tg) vs. observed frequency, v ohs for f2000 (solid, blue), f2000C (short- 
dashed, red), f250 (long-dashed, green) and f250C (dotted, magenta). We 
also show the corresponding 5Tb ((n 2 ) 1//2 /n), i.e. the temperature fluctu- 
ations if they were following the fluctuations of the density, normalized to 
their respective means (thin lines, same line types and colours for the cases), 
(bottom) the density (thin, solid, black) and differential brightness temper- 
ature fluctuations relative to their respective means (thick lines, same line 
types and colours for f2000, f2000C, f250 and f250C as in the top panel) 
vs. frequency. 



more than an order of magnitude larger than the relative fluctua- 
tions of the underlying gas density field. 

Although less of a requirement than in the case of the global 
step, a sudden decrease in the rms fluctuations should be easier to 
detect observationally, since other, possibly confusing sources of 
fluctuations, e.g. due to the instrumental setup or the foregrounds, 
are expected to change only slowly with frequency. Our f2000 sim- 
ulation does indeed shows a fairly steep drop over 15 MHz after the 
peak (but located in the FM band), whereas the other cases display 
a more gradual change over 20-30 MHz. 



7.2 Power spectra 

The quantity closest to the actual radio array observational data 
is the 2D angular power spectrum. It shows how the power of 
the fluctuations is distributed over the different angular scales. 
We follow the convention used in CMB analysis and construct 
the angular power spectrum expanded in spherical harmonics, 
[£(£ + 1)G/2tt] 1/2 , where £ = 2tt/0, with the angle in radi- 
ans. The results are shown in Figure 12 for simulations f2000 and 
f250 (solid lines). For reference we also show the corresponding 
angular power spectra of the underlying density field (dotted lines). 
These power spectra were constructed at full resolution without any 
beam- or bandwidth smearing, since these are instrument- specific. 
The smallest value of £ at which we calculate the power spectra 
corresponds to half of our box size, since below that there is an un- 
physical damping due to the finite simulation box size. The largest 
£ value shown corresponds to one computational cell and thus is 
related to our spatial resolution. 

The power spectrum at redshift z — 20.3 is identical to that of 
the underlying density field, since at that time there are only a few, 
small H II regions. In both reionization histories, evolution of the 
power spectrum proceeds by a rise at small scales (large £), shift- 
ing slowly to larger scales as reionization reaches the 50% point, 
and followed by a decline of the power as most of the universe gets 
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Figure 13. Angular power spectra [1(1 + 1)Q/2tt] 1 / 2 at z = 13.62 for 
f2000_406 and using bandwidths (top to bottom) 1 cell, 0.1 MHz, 0.2 MHz, 
0.5 MHz, 1 MHz, and 2 MHz. This illustrates the effect on the power spec- 
trum of integration over a given bandwidth. 



ionized. The declining power spectrum keeps peaking at approxi- 
mately the same £ value as when it was at its maximum, but the 
distribution becomes quite broad. Comparing the f2000 and f250 
cases shows that the f250 case peaks at somewhat smaller scales 
(i « 5000, « 4') than the f2000 case (£ « 4000, « 5'). The 
f250 distribution is also somewhat broader with a slightly lower 
peak value. 

This evolution of the power spectrum reflects the evolution of 
the H II regions. Initially they are small, but dominate the fluctua- 
tions, resulting in more fluctuation power at high £. The maximum 
then moves to lower £ as the H II regions grow with time. As we 
showed in Paper I this growth levels off at a certain size (~ 10 Mpc 
for f2000) because once an ionized bubble grows larger than that 
size it merges with many other bubbles to form much larger ionized 
regions. This reflects a typical correlation length between source 
clusters, convolved with the source efficiency (number of ionizing 
photons produced per atom). The last is reflected in the f250 power 
spectra peak being at lower scales (higher £). These very large bub- 
bles which result from the mergers of multiple smaller ones are 
rare, however, and locally they partly retain the shapes and sizes of 
the original smaller bubbles which merged (see e.g. Figure 1), so 
the typical fluctuation size at the scales under consideration (sub- 
degree) is still dictated by the typical size (~ 10 Mpc) of bubbles 
resulting from local source clustering. 

The power spectra magnitude grows strongly during the ini- 
tial, very patchy phase of reionization, in all cases reaching a max- 
imum of rsj 10 mK around the time of 50% ionization, in complete 
agreement with the results from the rms fluctuations we discussed 
above in Sect. 7.1. This is substantially over the power contained 
in the density field, about a factor 2 for the f2000 case, and 1.5 
for the f250 case. At late times (z < 12.3 for f2000) the power 
spectra are significantly depressed due to the small neutral fraction 
remaining, with the peak value reaching only up to 1-2 mK. The 
broadness of the peak indicates that the power is contained in ever 
larger structures. 

We do not show how beam smoothing affects the power spec- 



tra as this can be derived easily by multiplying the power spectrum 
with that of the beam (Fig. 5), and for detailed observational pre- 
dictions this should be done with the "dirty beam" of the particu- 
lar array under consideration. Both the Gaussian and compensated 
Gaussian beam will introduce an effective cut-off at the FWHM an- 
gle. Considering the multipole values for the peaks which we find, 
this means that planned experiments may be just sufficient to cap- 
ture the expected peaks. Our results are qualitatively similar to the 
analytical estimates in Furlanetto et al. (2004b,a). These authors 
also find that during reionization a peak develops in the angular 
power spectrum, and as reionization progresses it shifts to lower L 
However, our peaks tend to be more pronounced, and do not reach 
as low ts as in the analytical estimates. This divergent behaviour 
reflects the differences in their prescription for the growth of bub- 
bles and our simulations. 

The bandwidth smoothing can in principle be derived by ap- 
plying the appropriate window function to the 3D power spectra. 
Qualitatively, bandwidth smoothing takes away power at the high £ 
values since it will average of small bubbles along the line of sight. 
To illustrate the effect we show in Fig 13 how bandwidth smoothing 
affects the power spectrum at z — 13.62 for the f2000_406 case. 
The bandwidth varies from about 30 kHz (one cell) up to 2 MHz. 
Without smoothing the peak is at £ ~ 4000, moving to larger scales 
for larger bandwidths. However, both the magnitude and the posi- 
tion of the peak barely change for Ais ^ 0.2 MHz, staying close 
to the maximum value of 10.3 mK. For the maximum bandwidth 
considered, 2 MHz, the peak value drops by almost a factor of 2, 
to 5.46 mK, and occurs at £ ~ 3000. As expected, even modest 
smoothing washes out most of the small-scale power, but there is 
also some decrease of power on the larger scales, since we sample 
larger scales at which the fluctuations are lower. 

7.3 Beyond the power spectra: non-Gaussianity of the 21-cm 
signal 

The power spectra discussed in the previous section provide impor- 
tant information about the magnitude of the 21-cm fluctuations at 
different scales. However, they fall well short of a full description 
of the statistics of the 21-cm fluctuations, and the corresponding 
implications for high-redshift structure formation and reionization. 
Distributions with identical power spectra could have completely 
different properties and statistics. In Paper I we showed quantita- 
tively for the first time that the PDFs of both the ionized gas frac- 
tion and the ionized gas mass are generally strongly non-Gaussian, 
even more so than the underlying density distribution, which starts 
Gaussian, but develops non-Gaussian features due to the forma- 
tion of non-linear structures. Here we extend our results further by 
studying the PDF of the 21-cm differential brightness temperature. 

We calculate the PDFs using the same method as in Paper I, in 
cubes of 5, 10 and 20 h' 1 Mpc size (approximately corresponding 
to angular resolutions of 10', 5' and 2/5 respectively) and at two 
different redshifts for simulation f250, the first at the 50% ioniza- 
tion (z = 11.9) and the second during the late stages of reioniza- 
tion (z = 10.8). For the other simulations the PDFs are similar at 
the same ionization levels (which correspond to different redshifts 
due the different reionization histories). The results are shown in 
Figure 14 (solid lines), along with the corresponding Gaussian dis- 
tributions with the same means and widths (dotted lines). These 
Gaussian distributions are the ones one would construct on the basis 
of the power spectrum analysis. The PDFs are expressed in terms 
of deviations from the mean differential brightness temperature at 
that redshift. 
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Figure 14. Non-Gaussianity of the 21 -cm signal: PDF distribution of the 21 -cm signal from simulation f250 for z = 11.9 (left) and z = 10.8 (right), shown 
both at logarithmic (top row) and linear (bottom row) scale. The PDF were derived for cubical regions of sizes 20 h~ 1 Mpc (black solid), 10 h~ 1 Mpc (red 
solid), and 5 h~ 1 Mpc (blue solid). Also indicated are the Gaussian distributions with the same mean values and standard deviations (dotted, corresponding 
colours). 



These plots show that the PDFs of the 21 -cm signal are highly 
non-Gaussian. Considering for instance the 5 h~ x Mpc case (blue 
lines) at z — 11.9, there is a 5-10 times larger probability for find- 
ing holes (corresponding to ionized regions) deeper than —10 mK 
from the average signal (which is 16 mK at that redshift) than ex- 
pected from the Gaussian analysis. For the same 5 Mpc scale the 
peaks of 10 mK above the mean are also more frequent than the 
Gaussian statistics predicts, but higher peaks (> 12 mK) are actu- 
ally more rare. The same trends are seen at larger scales (10 and 
20 hT 1 Mpc), although to somewhat lesser extent. 

At the later redshift of 10.8, at which time the mean signal is 
5 mK, there are very large over-abundances, by factors of up to 30, 
of holes of —6 to —2 mK. Even more interestingly, strong peaks of 
10-20 mK above the mean signal are up to 1-2 orders of magnitude 
more abundant than the corresponding Gaussian predictions, and 
the distribution is even more skewed than at the earlier time. This 
over- abundance in the number of high peaks is roughly indepen- 
dent of the spatial scale considered, but occurs at different temper- 
atures, from rsj 10 mK for 20 h~ x Mpc regions, up to 16-18 mK for 
5 h~ x Mpc regions. 

This shows that the even very late into reionization, when most 
of the universe is ionized (for simulation f250 the mass ionized 
fraction is 84% at z = 10.8), there is still considerable signal avail- 



able in isolated features. The results in Sect. 6.1 already hinted at 
this, and the non-Gaussian statistics puts numbers to it. The same 
behaviour is seen in all of our simulations, and should increase 
the chances of actually detecting the redshifted 21 -cm line, even 
if reionization was mostly complete rather early. 



We can compare our results to previously published reioniza- 
tion PDFs. Ciardi & Madau (2003) studied a 20 /h Mpc simulation 
volume and presented the PDF at the resolution of their pixel size, 
much below the resolution of planned observations and the sizes 
we are addressing here, making a comparison difficult. Furlanetto 
et al. (2004a) calculated PDFs from their semi-analytical excursion 
set model. The z — 13 curve in their Fig. 4 can be roughly com- 
pared to our 5 Mpc curves at z = 11.9 since both refer to a state 
of rsj 50% ionization seen at ~ 2/5. In order to facilitate the com- 
parison, Fig. 14 also shows our PDFs plotted on a linear scale. We 
see that our results differ substantially from theirs, not being Gaus- 
sian as in their outside-in and uniform reionization models, and not 
showing the sharp cut-off of their inside-out models. 
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8 CONCLUSIONS 

We constructed the evolution of the redshifted 21 -cm emission for 
several reionization histories resulting from large-scale radiative 
transfer simulations in a volume of (lOO/i -1 Mpc) 3 . This volume 
corresponds to ~ 1 degree on the sky and ~ 5 — 10 MHz in fre- 
quency, sufficient for capturing the scales relevant to reionization 
and its 21 -cm signatures. We studied the various aspects of these 
21 -cm signatures, from the properties and statistics of individual 
bright features, to sky radio maps, to statistically-extracted signals. 
Our results can be summarized as follows. 

We do not find a very sharp global step in the 21 -cm sig- 
nal resulting from the neutral gas being ionized and disappear- 
ing. Instead, the transitions we observed are rather gradual, with a 
~ 20 mK decrement of the signal over ~ 20 MHz. A similar results 
if obtained for the change in rms fluctuations with frequency. As- 
suming well-behaved foregrounds, such a transition is well within 
the achievable instrument sensitivity and thus could still be de- 
tectable, though not as easily as a sharper transition would have 
been. 

The large-scale geometry of reionization is well-resolved in 
21 -cm sky maps even after beam smearing. Such maps would show 
structures in the neutral IGM, as well as the shapes of the HII re- 
gions are the sources, providing a unique window on early struc- 
ture formation. We also showed that subtracting the average signal 
from a large frequency band (~5MHz) should be efficient way to 
remove an idealized smooth foreground while leaving the 21 cm 
fluctuations essentially intact. 

Analyzing high-resolution LOS spectra we find that redshift 
distortions are considerable, and hence will influence the shapes 
of structures along the LOS, and should be taken into account 
when constructing 3D (position/frequency) observables from sim- 
ulations. There is a large variation between different LOS in terms 
of the structures encountered and the time when the epoch of over- 
lap is reached. If the foreground emission from galactic and extra- 
galactic sources can be effectively removed, the fluctuations seen 
along the LOS would provide detailed information about the un- 
derlying neutral density field. However, the most easily detectable 
objects would be local ionized bubbles, which are bordered by 
sharp transitions, where the 21 cm signal changes by tens of mK 
over very small intervals in frequency. We also presented for the 
first time simulated LOS spectra filtered with a realistic beam and 
bandwidth. These show that although small details are smoothed 
over and the signal variations are smaller than in the full-resolution 
spectra, the main structures remain clearly visible and should be 
detectable. 

Statistical measures, such as fluctuations and the two- 
dimensional power spectra show that the fluctuations always reach 
a maximum when the gas is ~ 50% ionized by mass, and that the 
power is boosted considerably (by a factor ~ 2) compared to the 
total density. At that maximum the power spectra show a clear peak 
at an angular scale of £ ~ 3000 — 6000, depending on the partic- 
ular reionization history. Earlier and later, the peaks are broader, 
indicating that more scales, both larger and smaller are involved. 
In order to be able to observe this peak the planned experiments 
should have a resolution of ~ 3' or better. 

The derived PDFs show that the 21 cm signal is strongly non- 
Gaussian at all scales. In particular, at late times the chance of find- 
ing individual bright features can be 10 or more times above the 
Gaussian prediction. At all times the brightest point in our simu- 
lation volume is 20-60 mK for beam/bandwidth smoothing typical 
for the planned LOFAR/PAST/GMRT-type observations. This is in- 



sensitive to the particular reionization history. It is therefore likely 
that even in the case of early reionization, individual bright peaks 
can be observed at frequencies above 120 MHz. 

Finally, we note that while we consider several specific reion- 
ization scenarios here, many of our conclusions appear to be quite 
generic, e.g. the significant boost of the 21 -cm emission fluctua- 
tions due to reionization patchiness, the fluctuations reaching maxi- 
mum at rsj 50% ionization by mass and declining thereafter, and the 
relative insensitivity of the magnitude and statistics of the brightest 
spots in the box. On the other hand, certain important features vary 
significantly among scenarios - e.g. at which redshift the fluctua- 
tions peak is reached. These depend on a number of still highly- 
uncertain reionization parameters. In the models discussed here we 
vary the most important of these parameters, the ionizing photon 
production efficiency of the sources and the level of gas clumping 
at small scales. As we discussed above, the ability to detect EOR 
signatures with any particular instrument would depend strongly on 
when these features occur. 
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